The Critical Properties of Two-dimensional Oscillator Arrays 
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We present a renormalization group study of two dimensional arrays of oscillators, with dissipa- 
tive, short range interactions. We consider the case of non-identical oscillators, with distributed 
intrinsic frequencies within the array and study the steady-state properties of the system. In 
two dimensions no macroscopic mutual entrainment is found but, for identical oscillators, critical 
behavior of the Berezinskii-Kosterlitz-Thouless type is shown to be present. We then discuss the 
stability of (BKT) order in the physical case of distributed quenched random frequencies. In order 
to do that, we show how the steady-state dynamical properties of the two dimensional array of 
non-identical oscillators are related to the equilibrium properties of the XY model with quenched 
randomness, that has been already studied in the past. We propose a novel set of recursion 
relations to study this system within the Migdal Kadanoff renormalization group scheme, by 
mean of the discrete clock-state formulation. We compute the phase diagram in the presence of 
random dissipative coupling, at finite values of the clock state parameter. Possible experimen- 
tal applications in two dimensional arrays of microelectromechanical oscillators are briefly suggested. 

PACS numbers: 64.60A-64.60ae,64.60.Bd, 64.60. Ht 



I. INTRODUCTION 

The study of the dynamical properties of large ar- 
rays of self-sustained oscillators with distributed intrin- 
sic frequencies is an interesting problem bridging non- 
equilibrium statistical physics and non-linear science. 
Large populations of interacting, self-sustained oscilla- 
tors are known to model a great variety of biological 
systems and several progresses have been made in the 
last decades [H, [3, If dense arrays of coupled self- 
oscillators, close to an oscillatory instability of the Hopf 
type, are expected to present interesting critical behav- 
ior, the analysis of such systems in the low dimensional 
case, i.e. when the range of the interactions is short, has 
been quite a challenging one and very little attempts have 
been made in this direction. The mean-field Kuramoto 
model of phase oscillators represents a very special case 
where one can predict in a clear way the macroscopic 
mutual entrainment (MME) properties of the system, but 
attempts to generalize such results to the low dimensional 
case are usually difficult and the majority of the studies 
of oscillator arrays have been devoted to the all to all 
case, and to networks with relatively high degree of con- 
nectivity. Meanwhile, the case of short range interactions 
is a very interesting one and have a lot of potential appli- 
cations, in particular considering the case of dissipatively 
and/or reactively coupled arrays of mechanical oscillators 
at micro and nano scales. 

The oscillator lattice problem, namely the finite- 
dimensional, nearest neighbour version of the Kuramoto 
model, has been attacked in the past by several authors 
[3, The main conclusions are that the synchroniza- 
tion properties one normally finds in the mean-field case 
are drastically changed by the short range nature of the 
connectivity, and that no synchronization is expected in 



two dimensions [1, H, @] ■ An interesting analysis based 
on real space renormalization group theory Q suggests 
an explicit expression for the lower critical dimension, 
below which one do not expect macroscopic mutual en- 
trainment, the latter being the dynamical analogue of fer- 
romagnetic order. A related situation occurs within the 
equilibrium statistical physics of spin systems with con- 
tinuous degree of freedom of the XY type, as already sug- 
gested in In the statistical mechanics of continuous 
spin models of the mean-field type the system has a spon- 
taneous magnetization and the out-of-equilibrium behav- 
ior is a very interesting and open problem These 
models relate closely to the Kuramoto model, where the 
additional difficulty of the intrinsic, quenched frequency 
distribution has to be taken into account, and where the 
system does not reach equilibrium but is self-driven in an 
out-of-equilibrium steady state regime. 

Both the dynamical and equilibrium properties of 
continuous spin models change importantly when short 
ranged, nearest neighbour interactions are considered. 
Rather than ordinary symmetry breaking, in the two- 
dimensional XY model algebraic order (AO) occurs. 
Namely the system has infinite correlation length, but 
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no magnetization at finite temperature 

When considering the case of low dimensional arrays of 
coupled oscillators, the problem is further complicated by 
the fact that the system is driven out-of-equilibrium, so 
that the standard tools of equilibrium statistical mechan- 
ics cannot apparently be exploited. On the other side, if 
the analogy between arrays of identical oscillators and the 
statistical physics of XY type of models with continuous 
degrees of freedom is a well known fact and attempts to 
relate the dynamical properties of low-dimensional arrays 
of oscillators to the classical theory of dynamical critical 
phenomena [l^, [l^l have been made, it is interesting to 
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see if the great deal of results obtained in the past for the 
finite dimensional XY model in the presence of disorder, 
could possibly relate to the problem of arrays of dissi- 
patively coupled, short ranged oscillators with a natural 
intrinsic frequency spread. We refer to this last prob- 
lem as the finite dimensional Kuramoto model or lattice 
oscillator problem Obviously the simplest case one 
may consider is the strictly two dimensional case, where 
the XY type of models are widely studied and where real 
space renormalization group theory has been capable of 
predicting a great deal of results in the past 

We will show that the formal analogy between identi- 
cal, dissipatively coupled arrays of oscillators and model 
A of non-equilibrium critical phenomena can be 
extended, in the two dimensional case, to the case of 
non-identical oscillators, so that the steady state non- 
equilibrium properties in the system may be studied as 
a function of the quenched intrinsic frequency distribu- 
tion. This will allow us to describe the steady state out 
of equilibrium properties of two dimensional oscillator ar- 
rays by mean of classical statistical mechanics. The pos- 
sibility to consider non-identical oscillators is essential, 
in that a finite width of the intrinsic frequency distri- 
bution is the fundamental parameter one consider in all 
Kuramoto type of models. We will also include a physi- 
cal temperature corresponding to white noise in the array 
and consider the case of bimodal frequency distributions, 
options also commonly considered in the context of the 
mean- field Kuramoto model f§\ . 

In this respect we will show that the steady state prop- 
erties of a two dimensional array of dissipatively coupled 
non-identical oscillators can be related to the equilibrium 
properties of an effective XY model in the presence of 
random Dzyaloshinskii-Moriya (DM) interactions. The 
random (DM) interaction effectively induces phase cant- 
ing between neighbouring oscillators, and can be related, 
as we will show below, to the quenched distribution of 
the native intrinsic frequencies. It is then argued that 
the steady state non-equilibrium properties of two di- 
mensional arrays of dissipatively coupled non-identical 
oscillators (in the phase approximation) are related to 
the equilibrium properties of the two dimensional XY 
model in the presence of disorder, where again a great 
deal of real space renormalization group theory results 
are available, even though the effect of quenched ran- 
dom interactions on the classical properties of the two 
dimensional XY model is still a matter of debate within 
the condensed matter physics community. A similar ap- 
proach to describe the non-equilibrium properties of two 
dimensional superconducting arrays with external cur- 
rents by mean of the equilibrium properties of an effective 
XY model has been considered in the past [Tsj . 

Before we begin our analysis, it is important to stress 
that the properties of the Kuramoto model in low dimen- 
sions have already been discussed in the past 0, H, 0, 0] ■ 
In particular one could argue at this point that in two 
dimensions no macroscopic mutual entrainment (MME) 
is expected, so that collective synchronization is simply 



not present. This argument reminds a similar one about 
the very nature of the critical properties of the two di- 
mensional XY model [l^, [HI . If magnetization is known 
not to be expected at finite temperatures, the two di- 
mensional XY model is known however to present crit- 
ical properties of a special type (AO). We will show in 
the following that the dynamical analog of algebraic or- 
der is present in dissipatively coupled two dimensional 
arrays of oscillators (and we might refer to it as alge- 
braic synchronization). Dissipatively coupled oscillators 
in two dimensions are characterized by a phase transi- 
tion of the Berezinskii-Kosterlitz-Thouless (BKT) type, 
as we will show in this article, even for non-identical os- 
cillators, namely when one consider the non trivial case 
of distributed intrinsic frequencies. The main question is 
how algebraic order (AO) reacts to the disruptive effect 
of the intrinsic frequency distribution, and if a steady 
state algebraic synchronization regime exists in two di- 
mensional arrays of oscillators. 

We will address these questions by evaluating the 
renormalization group recursion relations, within the po- 
sition space Migdal-Kadanoff (MK) approximation, for 
the case of an XY model with quenched random inter- 
actions induced by effective Dzyaloshinskii-Moriya (DM) 
interactions. We consider the standard formulation of 
the XY model within the (MK) approach [13], adopt- 
ing the discretization scheme introduced in reference |18| 
and proposing a novel set of recursion relations that will 
allow us to include both random exchange and random 
(DM) interactions, consistently with the usual methods 
of real space (MK) renormalization group of disordered 
spin systems [l^. We will be considering a renormal- 
ization group rescaling length 6 = 3, to treat exchange 
and random (DM) interactions of random opposite signs 
on equal footing [20] , effectively generalizing the position 
space renormalization group (PSRG) recursion relations 
of the two dimensional XY model to the case of random 
exchange and (DM) interactions. Together with the novel 
set of (PSRG) recursion relations, we present prelimi- 
nary results obtained by means of an algorithmic method 
capable to implement the recursion relations and evalu- 
ate the corresponding phase diagram, for different model 
systems, corresponding to different forms of disorder, by 
fixing the proper choice of the initial conditions in the 
renormalization group transformation. Choices we will 
consider ranges from the XY spin glass (XYSG), where 
randomness is in the exchange interaction and no (DM) 
interactions are present, the two dimensional ferromagnet 
with Dzyaloshinskii-Moriya interactions (XYDM) prob- 
lem, where randomness is in the (DM) term and no ran- 
dom exchange interaction term is present, and ultimately 
the random gauge glass (RGG) problem, as an impor- 
tant general case where both type of randomness are 
present simultaneously, so that gauge invariance is re- 
stored. Our MK renormalization group approach will be 
general enough to describe any of the models above, by 
a proper choice of the initial conditions in the renormal- 
ization group flows. 
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The proper initial condition of primary interest, i.e. 
the one we wish to consider to be able to predict the 
type of behavior one expects in two dimensional arrays 
of dissipatively coupled arrays of non-identical oscillators, 
is, as we now show, of the (XYDM) type, even though 
randomness is the dissipative/exchange coupling term is 
also relevant when discussing two dimensional oscillator 
arrays. 

Our results do not restrain to the oscillator array prob- 
lem only, but are general enough and pertain to differ- 
ent type of systems that have been considered in the 
past within the condensed matter theory community, and 
where it is relatively common to consider the properties 
of the XY model in the presence of quenched random 
interactions of various types. 

The shape of the phase diagram of the two dimensional 
XY model in the presence of quenched random (DM) 
impurities is the subject of a long standing debate [2l|, 
[22 . 23, 24] we here address and show to relate also to the 
study of dissipatively coupled, two dimensional oscillator 
arrays. 



II. THE MODEL 

The dynamical behavior of oscillator arrays in the 
vicinity of a supercritical Hopf bifurcation can be de- 
scribed in terms of complex amplitude equations, related 
to the amplitude and phase of each oscillator. In the 
general case of reactive and dissipative coupling, as well 
as considering the presence of stiffening non-linearities 
of the Dufhng type, and in the case of an intrinsic fre- 
quency distribution spread within the oscillators, one can 
write the following equations for the complex amplitude 
[25', 'Sg'I, where, as mentioned above, we add a finite tem- 
perature due to brownian fluctuations within the array 
so that {■qk{t)-nk'{t')) = 2TSk,k'S{t - t'). We restrict the 
range of the interactions to neighbouring oscillators on 
the square lattice. 

dAk 

—— = {^j. + iLUk)Ak - + ia)\Ak\'^Ak 
dt 

+ {J + iR) Y,{A,- Ak)+'nk (1) 

We will assume that the oscillators are not identical 
so that quenched intrinsic frequencies u>k are considered, 
with a given (e.g. bimodal) distribution P{Ld). The two 
parameters and 7 are related to the self-sustaining 
mechanism that drives each oscillator; for oscillators of 
the Van der Pol type we can choose /i = 7. The coupling 
J is the dissipative coupling and we here assume to be 
uniform and short ranged. In what follows we will also 
consider the possibility of a random dissipative/exchange 
coupling. Finally the reactive term i?, which is related 
to a mechanical coupling between the oscillators, is also 
short ranged and eventually distributed around an aver- 
age value Ro- We note at this point that the interesting 



case of reacti yely coupled oscillators has been considered 
in reference [27|, in the mean- field approximation. We 
neglect this last term, and rewrite the above equation as 

^ - (A^ - irl)ru + J 5^ r, (1 - cos(0, - Bu)) + vl 
jeCk 

^ = [iVk - akrl) + JJ2 ^j/^k sm{e, - 9k) + 4 (2) 

Including a reactive term means one should include 
RJ2j£Ck sin(0j — 9k) to the right end side of the first 
equation and RJ2j<£Ck ^j/'''k{cos{9j — 9k) — 1) to the right 
end side of the second. We consider in what follows the 
isochronous case, assuming that the Duffing parameter 
ak can be neglected, as in the original formulation of [3]. 
The magnitude of the complex amplitude crosses over to 
one as soon as the value of /i, related to the Van der 
Pol strength of the autonomous oscillators, is properly 
tuned. Assuming the width of the intrinsic frequency 
distribution to be relatively high peaked around the av- 
erage value ujo, the above equation reduces to the phase 
approximation, and one writes 

^=c., + j5]sin(0,-0fe) + r;^ (3) 

where Ck are the neighbouring sites to k. In the sim- 
ple case of identical oscillators one recover the Langevin 
equations for the XY model. Differently, we are back to 
the oscillator lattice problem. In other terms, as men- 
tioned in the introduction, the dynamical behavior of 
identical oscillators, reduces, up to a global redefinition 
of the complex amplitudes, to the coarsening dynamics of 
the XY model, so that the above equation ([T]) reduces to 
model A of critical dynamical phenomena [UH^, and 
the formal analogy between a critical Hopf bifurcation 
and phase transition theory can easily be drawn in this 
case. The dynamical behavior of non-identical oscillators 
is certainly harder to solve, as usually happens when deal- 
ing with models in the presence of quenched randomness. 
We are effectively studying the low-dimensional version 
of the Kuramoto model, and the first question is how the 
random frequency term affects the onset of algebraic or- 
der (AO) we know to be present in the case of identical 
oscillators. We will address these questions in the rest of 
the paper, showing how the non-equilibrium steady state 
properties of the problem can be related to the equilib- 
rium properties of an effective XY model with random 
quenched interactions. Before doing so, it is interesting 
to review how (and if) numerical simulations might be 
useful to answer, at least in part, this question before we 
return to a real space, renormalization group calculation 
in the rest of the paper. We clearly expect that disorder, 
in combination to the finite dimensional character of the 
connectivity, might result in numerical simulations to be 
converging very slowly, a situation we faced and we will 
discuss in the following section. 
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FIG. 1: Helicity Modulus computed for system sizes L — 
5, 10, 1, 20, 25, 50 in the two dimensional XY model without 
disorder, by mean of direct Runge-Kutta integration meth- 
ods. The Helicity modulus jump becomes steeper for increas- 
ing system sizes. We could extrapolate a critical transition 
value of Tkt — 0.89, in agreement with previous, accurate 
estimates [331. 



III. NUMERICAL SIMULATIONS 

In the current assumption of dissipative coupling and 
within the phase approximation we discussed above, 
we shown that the amplitude equations for the two- 
dimensional array reduce to the equations of motion of 
the XY model, with the inclusion of a quenched ran- 
dom frequency term. How are the critical properties of 
the XY model affected by the random frequency term? 
We firstly answer this question numerically, by means 
of direct integration methods. We considered a gaussian 
form for the intrinsic frequency distribution and solve the 
Langevin dynamics explicitly, via Runge-Kutta methods 
for a system of size N = L x L, following the standard 
methods discussed in the literature [13, ISH. We also 
checked that the same algorithm, in the case of mean- 
field interactions and when random frequency terms are 
included, reproduces the results of the Kuramoto model. 
We also check that we were able to obtain the expected 
results in the case of identical oscillators in two dimen- 
sions. Even though the system does not order at finite 
temperature, phase correlation effects are present and are 
quantified by finite values of the helicity modulus [s^ ]. 
as predicted by the Kosterlitz Thouless renormalization 
group theory. The numerical results we here show are 
consistent with the expected (BKT) transition [ll|, |3^ . 
The algebraic decay of the correlation function can be 
measured and the weak violation of universality critical 
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FIG. 2: Critical exponent 77. We computed explicitly the 
correlation function at increasing values of the temperature 
T < TkT, observing as expected, algebraic decay. The es- 
timated value of the exponent is consistent with the value 
77 = 1/4, expected at T = Tkt [H. 



phenomena is seen, as expected, with the critical expo- 
nent being close to the expected value 77 = 1/4 at the 
transition point (see Fig. [2]). In the presence of disor- 
der, as usually occurs in low dimensional systems, things 
become difficult from the perspective of numerics. We 
clearly observed that for each given quenched realization 
of the intrinsic frequency distribution, the system reaches 
a steady state, i.e., for any realization of the disorder, we 
observed that the Runge-Kutta algorithm reaches a time 
independent state. We computed the average helicity 
modulus over several realization of disorder, for a given 
system size ranging, as in the pure case, from L = 5 to 
L = 50, for two dimensional arrays of size L x L. We 
observed however that the results we obtained from nu- 
merics in the presence of disorder are not that clear, the 
sample-to-sample fluctuations becoming very strong as 
the width of the intrinsic frequency distribution grows. 
For very small system sizes we observe the helicity mod- 
ulus to be finite below the critical point, that decreases 
with the disorder strength, as in Fig. [31 However for 
systems of size as small as N ^ L x L = 16^ it becomes 
very difficult to have sample-to-sample fluctuations un- 
der control, so no conclusive results were obtained for 
these system sizes, even after a run of several weeks on 
a standard workstation. Differently, in the mean field 
case, where each oscillator is coupled to all others with 
equal strength, we were able to control sample-to-sample 
fluctuations and an average value of the Kuramoto order 
parameter was easily obtained, consistently with theo- 
retical predictions, up to relatively large system sizes 
N = 10^. To conclude this section on the direct in- 
tegration approach of ([3]), we want to comment further 
the results we obtained in the presence of disorder. Our 
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FIG. 3: Helicity Modulus computed for system size A'^ — 
16 X 16. For this small system size we were able to con- 
trol sample-to-sample fluctuations and estimate the helic- 
ity modulus for increasing values of the disorder strength, 
A = 0.05,0.1,0.2,0.3, deflned as the width of the gaussian 
distribution of the intrinsic frequency distribution P{i^)- 



findings seems to indicate that the location of the (BKT) 
transition decreases for increasing values of the disorder 
strength, as in the inset. On the other side, in the low 
temperature region and for system sizes L > 32, we do 
not see the helicity modulus to reach a constant finite 
value as for the small system L = 16 analyzed in Fig. 
[3l In other terms, if algebraic order is seen to occur at 
finite temperatures, for finite, small values of the intrin- 
sic frequency width distribution, numerics cannot state 
in a clear manner whether such transition extends all the 
way down to zero temperatures. The helicity modulus 
sample-to-sample fluctuations diverge in the low temper- 
ature region, so unique conclusions on the form of the 
phase diagram, in the space of temperature and disorder 
strength, cannot be reached by mean of direct integra- 
tion methods. We believe that the difficulties we just 
discussed do not really depend on the specific method 
we used, but simply, as often occurs in low dimensional 
systems, the presence of disorder makes numerical simu- 
lations a very hard task. 

In this section we have seen how direct integration 
methods cannot state in a clear manner how algebraic 
order is affected by the presence of random quenched in- 
trinsic frequencies. For this reason we will address the 
same question within real space renormalization group 
theory in what follows. 



IV. REAL SPACE RENORMALIZATION 
GROUP THEORY 



The effect of quenched random interactions on the 
Kosterlitz-Thouless type transition is an old and inter- 
esting problem that applies to a variety of different phys- 
ical problems. If one normally expects that infinitesi- 
mal bond randomness do not affect the transition, it has 
not been clarified in a conclusive way how other types 
of quenched random interactions, namely site disorder 
and/or random interactions of the (DM) type affects al- 
gebraic order. Real space renormalization group consid- 
erations suggests that at finite temperature the (BKT) 
transition is not destroyed [2l|. However the effect of 
quenched random (DM) impurities on the low temper- 
ature behavior of the system is usually expected to de- 
stroy (BKT) order, and a reentrant phase diagram has 
been predicted long ago, [2l|, 111, [s^l . On the other side, 
the validity of such results have been questioned in the 
recent past [23), i38|, and a Migdal-Kadanoff set of re- 
cursion relations have been considered |l3 , possibly sug- 
gesting that algebraic order is stable against disorder at 
finite, small temperatures so that no reentrance would 
be present. Before returning to this point, we want to 
stress why this question is also relevant when studying 
the effect of a random intrinsic frequency spread in the 
two dimensional array of dissipatively coupled oscillators 
we consider in this work. In the first part of this section 
we show that disorder in the intrinsic frequency distri- 
bution of the two dimensional oscillator array problem 
is related to the classical types of quenched random in- 
teractions that have been introduced in the past in the 
context of the XY model. Once this point is clarified and 
once we will explain how the non-equilibrium steady state 
properties of two dimensional arrays of dissipatively cou- 
pled array of oscillators can be understood studying the 
equilibrium properties of an effective XY model in the 
presence of quenched random interactions of the (DM) 
type, we will critically reconsider the Migdal Kadanoff 
position space renormalization group calculation of ref- 
erence [23] to establish/check whether and how random 
(DM) interactions affects algebraic order, proposing a 
new set of recursion relations to evaluate the correspond- 
ing phase diagram. We hence suggest the correct recur- 
sions that can answer the question of how quenched in- 
teractions affects (BKT) order according to the (MK) 
approximation, a method that has been shown to de- 
scribe low dimensional systems with quenched disorder 
in a simple and effective way. We believe that the phase 
diagram of the XY model, in the presence of random 
(DM) interactions within the MK approximation is con- 
sistent with the reentrant behavior originally predicted 
within real space renormalization group theory [2l[. We 
suggest however that reentrance reduces at very low tem- 
peratures, so that the phase boundary intercept the zero 
temperature axis at finite values of the disorder strength 
(with vertical slope). In order to obtain the phase dia- 
gram, we write generalized recursion relations for a dis- 
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cretized clock state model, following the work of reference 
[l^ . including bond and (DM) randomness in a consis- 
tent way, at any finite Q values of the clock state param- 
eter. Our recursion relations reduce to the well known 
random bond Ising model for values of the clock state 
parameter Q = 2 and reduces to the (MK) recursion 
relations of the pure XY model at large values of the 
clock state parameter Q, when no disorder is present. 
Without loss of continuity we are able to interpolate be- 
tween the Random Bond Ising model (RBIM) [s^, the 
pure XY model at large values of Q and no disorder as 
well as the general case of both exchange and (DM) ran- 
dom interactions being present, ultimately proposing a 
novel framework to study the random gauge glass model, 
where both exchange and (DM) interactions coexists and 
are related by the proper initial condition in the renor- 
malization group flows. 

The relevance of these models and of the above ques- 
tions to the oscillator array problem can be explained as 
follows. Let us consider the following effective hamilto- 
nian 

-pn = Y, J^S.Sj + J2 Drjz ■ X (4) 

{i,3) {i,3) 

In the classical formulation of the two dimensional XY 
model one expands around small values of subsequent 
phase shifts and then imposes the periodicity of the phase 
variables, leading to the canonical Coulomb gas formu- 
lation [40|, so the Kosterlitz Thouless renormalization 
group equations can be derived. Similarly, expanding the 
second term related to the random (DM) interaction, an 
effective Coulomb gas with random dipole moments can 
be obtained [2lj. Following these steps and averaging 
over the disorder either with the replica method or al- 
ternatively [21I, [3^ recursion relations of the Kosterlitz 
Thouless type that include the presence of disorder have 
been obtained and a reentrant phase diagram in the space 
of temperature and random strength has been predicted. 

The above hamiltonian, considering the case of no dis- 
order being present in the exchange interaction Jij = Jo 
in eq. usually referred as the two dimensional fer- 
romagnet with Dzyaloshinskii-Moriya (XYDM) interac- 
tions, is characterized by equation of motions that, to 
leading order in the adjacent phase shift 0i — 0j corre- 
sponds the the two dimensional Kuramoto model, as a 
simple calculation of the equations of motion reveals (see 
e.g. eq. (9.7) in reference [57]). We note that the in- 
trinsic frequencies uji are related to the quenched dipole 
distribution according to Ui = X]je£ lij^ where qij are 
the effective dipoles related to the original (DM) interac- 
tion. 

The system of oscillators with distributed, quenched 
frequencies has been related to the case of identical os- 
cillators with an effective quenched dipole field, that in- 
duces a relative frequency mismatch. According to the 
above hamiltonian ^ , and assuming no disorder present 
in the exchange interaction, the problem relates, as one 



can check explicitly writing the Langevin equations as- 
sociated with (HI, to the lattice oscillator problem in two 
dimensions. 

It is then reasonable to ask how to write the proper re- 
cursion relations within the the Migdal Kadanoff approx- 
imation, corresponding to the above hamiltonian system 

O- 

The equilibrium properties of model ^ relates to the 
steady state solution of the two dimensional oscillator 
model we originally intended to study. In the next sec- 
tion we present a study of the two dimensional XY 
model with quenched exchange and/or random (DM) in- 
teractions within the MK approximation. Not only our 
finding will be relevant to the oscillator array problem 
but to many other physical problems to which the two- 
dimensional XY model in the presence of disorder of the 
type dictated by equation (j4]) is known to relate. 

V. THE MIGDAL KADANOFF METHOD 

In the classical formulation of the pure XY model 
with the Migdal Kadanoff renormalization group scheme 
[l^ . renormalization group recursion relations for the 
effective coupling of combined bonds, within a Fourier 
mode formulation are considered. At low temperatures 
the approximation recovers effective algebraic order, the 
potential flows to a Villain form [13] and the system 
is characterized by an infinite correlation length, even 
though has vanishing magnetization. The only drawback 
of the (MK) approximation is that after a sufficiently 
large number of renormalization group steps, the recur- 
sion relations flows to a high temperature disordered sink 
at any finite temperature so that true fixed line behav- 
ior is not present in the strict sense. The number of 
iterations required for this to occur is however diverg- 
ing below the transition point so that effective algebraic 
order is successfully described by the approximation. It 
is then reasonable to ask how to reformulate the classi- 
cal (MK) approach to the XY model in the presence of 
quenched randomness, due to the overall success of the 
(MK) approximation in the context of disordered Ising 
type of models. 



A. The discrete Clock State 

Meanwhile a second, interesting formulation of the 
(MK) approximation for the two dimensional pure XY 
model has been suggested in the past 18]. Rather than 
the usual Fourier version of the (MK) recursion relations, 
one considers a discrete clock state model, at integer fi- 
nite values of the clock state parameter Q. The above 
formulation has been given for the pure XY model in two 
dimensions. Attempts in order to include the presence of 
(DM) interactions have been made [H, |4l|, |4j but we are 
not aware of any study where correct recursion relations 
have been written for this case. Another delicate point is 
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how one treats the recursion relations in the case of dis- 
order (both in the exchange and (DM) interactions). We 
treat disorder in the same way it has been considered in 
the past for the (RBIM), namely within a binning proce- 
dm'e, avoiding to resort to random pool methods |38l.l4^. 
our method being simply a deterministic one (despite we 
are dealing with disorder). Our recursion relations and 
algorithmic method reduce to the calculation we did in 
the past for the (RBIM) for values of the clock state pa- 
rameter Q = 2. The computational effort grows linearly 
with Q. As we will see for the simple case of the pure 
model, effective algebraic order is observed already for 
values of the clock state parameter as small as Q = 32. 
In order to treat random exchange and (DM) interac- 
tions of both sign on equal footing we need to consider a 
renormalization group rescaling length 6 = 3. This last 
choice makes the algorithmic solution of the recursion re- 
lations in the presence of disorder quite slow, since triple 
convolutions have to be considered. A similar issue has 
been discussed already for the simple case of the Q = 2 
(RBIM) li]. 



B. The Recursion Relations 

As mentioned above we consider the case of discrete 
phase angles 6i = 2'Kqi/Q, = 0, • • ■ , Q — 1. 
The hamiltonian reads 

- /3H = V Jy cos{e, -e,) + y] D,j sin(0, -Oj). (5) 



Under renormalization group transformation, we con- 
sider two generalized couplings between neighbouring 
sites Jij and Dij independently. The renormalization 
group recursion relations read 



D'{q) 

G'iQ) 



-\og{R{Q,q)R\Q,q))-G'{Q) 
\og{R{Q,q)/R\Q,q)) 

— ^log(i?(Q,g)i?t(Q,5)) 



(6) 



where interaction terms are always considered modulo 
Q, and where the captive renormalization group constant 
imposes the constrain J' {Q, q) = and a similar con- 
strain for D'{Q,q) equally applies. The renormalization 
group polynomials, for the case of a rescaling length fac- 
tor b — 2 read: 



Q-i 



RiQ,q)=Y.^''"' 



{l) + J23(q+l) + Di2{l)-D23{q+l) 



1=0 

Q-1 



RHQ,q)^Y. 



Jl2{l) + J23(q+l)-Di2il)+D23{q+l) 



where the "bond-moved" exchange interactions are 



71=1 



(8) 



(7) 



1=0 



together with a similar expression for the (DM) interac- 
tions Dij . These recursion relations might seem at first 
sight rather similar to the ones discussed in the past [l^ , 
however a important symmetry property is included in 
the above recursion relations, that was not discussed in 
previous formulations. A set of recursion relations of the 
(MK) type, where the two interactions were treated sep- 
arately as above, was discussed, but only in the so called 
harmonic approximation [4l|, |42, [43] . 

The above recursion relations will be firstly discussed 
for a renormalization group rescaling length 6 = 3, in the 
absence of disorder, so we can check that the onset of 
effective algebraic order described by the (MK) method 
is properly recovered within the discrete clock state for- 
mulation and we can determine the effect of an uniform 
(DM) interaction on the pure XY model. We will then 
consider the case of small clock state Q values, and in- 
clude the presence of disorder, checking that the above 
recursions reproduce, as expected, the phase diagram for 
the (RBIM) at Q = 2 precisely, before returning to the 
ultimate issue of large values of Q and disorder being 
present, corresponding to the XYSG, XYDM and RGG 
models discussed above. 

Attempts in this direction were already considered in 
[3, Hi except we are now able to incorporate the pres- 
ence of (DM) interactions in a consistent way. The re- 
cursion relations despite their simplicity, are the only 
possible ones that properly reflects the symmetries of the 
original hamiltonian, as a direct analysis reveals. Note 
that, as also explicitly stated in reference [l^, the sym- 
metry properties of the generalized potential there pro- 
posed were lost under renormalization group transfor- 
mation in the case of (DM) interactions being present, 
something that should simply not occur. The fact that 
the recursion relations written in the past were not com- 
plete is probably the reason behind the erratic behaviour 
of the renormalization group flows, also reported in 3^, 
whenever interactions of the (DM) type were considered. 
Any conclusion on the shape of phase diagram of the 
XYDM problem based on these recursion relations [2^ 
should then be taken with care. 

Differently, in our method the potential maintains its 
symmetry properties under renormalization group trans- 
formation (see e.g. Fig. [5]), as it should be under any 
(symmetry) group transformation. Once the novelty of 
the above recursions is understood, we need to discuss 
how to implement them algorithmically, in the presence 
of quenched random interactions. The way to treat dis- 
order from the algorithmic point of view is a delicate 
issue we should also discuss. As explained above we use 
a deterministic method based on a binning procedure, 
that is known to respect the symmetries of the renormal- 
ization group transformation, as already discussed in the 
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FIG. 4: The fixed Villain potential observed for the pure XY 
model, at temperature value T = 0.5 for a renormalization 
group rescaling length & = 3, in the absence of (DM) interac- 
tions. Crosses shows the fixed potential in the discrete scheme 
with Q = 32, while the continuous line shows the results ob- 
tained for the value of the clock state parameter Q = 512. 

context of the (RBIM) , that corresponds to Q = 2 in the 
above recursions. 



VI. IDENTICAL OSCILLATORS 

We first show the results we obtained for the pure 
model (and b = 3), at large finite values of Q, checking 
that the results of the (MK) approximation in the XY 
model are properly reproduced within the discretization 
scheme here considered and discussing the role of (DM) 
uniform interaction terms, before returning on the issue 
of disorder, we will be able to include via the proper 
choice of the initial conditions of the renormalization 
group flows. 

We do not report explicitly the form of the renormal- 
ization group polynomials for the case of 6 — 3, but 
they clearly include a double sum, rather than the sin- 
gle sum in eq.(U|j as usually discussed in the context 
of the (RBIM) [3^, as one needs to decimate an even 
number of sites at each KG step in order to treat interac- 
tions of both signs in the same manner, so to produce a 
phase diagram that will have the proper symmetry prop- 
erties. Going back to the pure case, when no randomness 
is present, we expect that the location of the Kosterlitz 
Thouless transition, within the (MK) approximation, de- 
pends on the decimation parameter b, as can be seen in 
Fig. [SI showing the number of iterations needed before 
the above recursions ([6]) flow to the high temperature 
disordered sink. Note that in the case of the pure XY 
model without (uniform) (DM) terms, one recovers the 
usual Villain potential 17] for any value of the clock state 



FIG. 5: The two components of the complex fixed poten- 
tial, observed for the pure XY model, at temperature value 
T = 0.5 for a renormalization group rescaling length 6 = 3, in 
the presence of uniform (DM) interactions, obtained for the 
clock state parameter value Q = 512. Clearly, under renor- 
malization group transformation, the exchange interaction is 
symmetric around n, while the (DM) generalized interaction 
is antisymmetric. Fix-line behavior as in the case where (DM) 
interactions are absent is observed. 



parameter Q > 16. We show, for the pure case, both the 
fixed Villain potential we obtain for Q = 32 (blue crosses) 
as well as Q = 512 (continuous red line), in Fig. [H 

These results shows that the clock state discrete ap- 
proximation is a very robust one, in that effective alge- 
braic order is found at relatively small values of the clock 
state parameter, and there is no need to increase the val- 
ues of Q of the order of 10"^, as done in reference [l8| . 
This last remark will be crucial when dealing with the 
disordered case, even though the computational effort of 
the algorithm we will introduce in the next section ul- 
timately scales in a linear way with Q, while it scales 
as a cubic power of the number of bins we will use to 
discretize the quenched probability distribution(s) of the 
exchange and (DM) interactions at each renormalization 
group step. 

In the presence of (DM) interactions the complex po- 
tential converges to a symmetric and an antisymmetric 
part ,as in the original interaction eq. Note that the 
symmetry (antisymmetry) of the two terms in the initial 
hamiltonian are preserved by the above recursion rela- 
tions ([6]), and that is because we wrote two distinct re- 
cursion relations for each term, since the renormalization 
group polynomials have to be regarded, in the presence 
of (DM) interactions, effectively as complex quantities. 
As in the case of the pure XY model, effective algebraic 
order is observed below the (BKT) transition point, that 
do not changes in the presence of uniform (DM) interac- 
tions, as expected. 

Before we conclude this section on the (MK) clock state 
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FIG. 6: The Number of Iterations for the pure XY model 
with 6 = 2 and 6 = 3. The results shows that effective al- 
gebraic order is observed in that the increase in the iteration 
step number is rather sharp around the effective critical tran- 
sition point. The results seems to suggest that the location of 
the transition point depend, as occurs in Ising type of mod- 
els within the (MK) approximation, on the renormalization 
group rescaling length. 

approach to the pure XY model with uniform (DM) in- 
teractions, we want to add that, as stressed above, the 
location of the (BKT) transition do not changes impor- 
tantly for values of the clock state parameter Q > 32, 
so we can safely consider the disordered case considering- 
finite values of Q in this range of values. 

For the case of the clock state parameter Q = 4 we ob- 
serve a ferromagnetic transition point at Tc ~ 2.078 that 
decreases for increasing values of the relative interaction 
strength D/J, as can be seen in Fig. [7l 



VII. THE ROLE OF DISORDER 

In order to test our recursion relations and algorith- 
mic methods, we begin reconsidering the simple (RBIM) 
case, that is included in the above equations at values 
of the clock state parameter Q = 2, 4, for vanishing val- 
ues of the (DM) interaction strength D. As expected 
we recover the reentrant phase diagram that has been 
discussed at length in the past. Clearly, a,t Q — 2,4, 
gaugein variance is respected, so we plot the Nishimori 
line [301 and note that, as already found in the past, the 
phase boundary reaches the highest value of the antifer- 
romagnetic bond concentration on this line. Note that, 
at _D = only, the Q — 2 and Q = 4 models coincide up 
to a redefinition of temperature so we discuss the latter 
case only. Interestingly, high precision measures at very 
low temperatures shows that reentrance, that is clearly 
present as reported already [39] , diminishes gradually at 
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FIG. 7: The location of the critical transition point for the 
clock state parameter value Q = 4, when no disorder is 
present, and where uniform (DM) interactions are considered. 
The a;-axis shows the relative strength of the two interactions, 
as chosen in the initial conditions or the recursions ((6)1. 



very low temperatures, so that our results are consistent 
with the expectation that the phase boundary is verti- 
cal, but at low temperatures only. Simply this do not 
occur all the way up to the Nishimori line, according to 
our findings. This interesting reentrant behaviour that 
gradually reduces to intercept the zero temperature axis 
with a vertical slope has also been observed recently in 
the context of (gauge invariant) Gallager codes, and we 
conjecture that the qualitative behaviour of the phase 
boundary we are here discussing applies to all spin sys- 
tems for higher values of Q, whenever gauge invariance 
is respected, except the ferromagnetic phase is replaced 
by a (BKT) phase for values of the clock state parameter 
Q > 16. It is important to mention at this point that 
any initial condition in the (RG) flows for Q = 2,4, close 
to the boundary and above the Nishimori line, will flow 
to the finite temperature unstable fixed point, while any 
initial condition below the Nishimori line will flow to the 
strong coupling low temperature fixed point already dis- 
cussed in the past [39| . a phenomenon known as strong 
violation of universality. 

The phase diagram above (Fig. [8]) clearly do not relate 
to the original XY model we planning to consider, but 
we performed this calculation to test the algorithm that 
implements the recursion relations in the presence of dis- 
order and we will be able to consider large values of Q of 
the (DM) type at a later stage, changing without loss of 
continuity the parameters in the initial condition of the 
RG flows. Moreover, it is quite interesting to study the 
clock state model at small finite values of the parameter 
Q with disorder, since it is a well defined problem that 
interpolates between the Ising and XY type of models. 
As found in the pure case, when disorder is present in the 
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FIG. 8: The Phase Digram for the clock state parameter 
value Q — 4, without (DM) interactions and in the presence 
of bimodal exchange interactions, with an antiferromagnetic 
bond concentration p. The inset shows the reentrant behavior 
of the phase boundary around the multicritical point. Gauge 
invariance is present in the model and so is reentrance. 



exchange interaction, we observe ferromagnetic order at 
low temperatures for the values Q — 2,4,8, while effec- 
tive algebraic order appears at values of the clock state 
parameter Q > 16. This means that the recursion rela- 
tions ^ flows to a ferromagnetic fixed point for values of 
Q < 8, while at Q > 16 we observe the usual quasi- fixed 
line behaviour discussed in the context of the XY model 
( see Fig. [TO)) , and we do not expect things to change im- 
portantly for higher values values of Q, as we have seen 
already for the pure case, meaning that the discretiza- 
tion scheme converges rapidly to the original problem of 
continuous degrees of freedom. For the case of Q = 8, 
and when DM interactions are not present we computed, 
for a bimodal exchange interaction, the phase diagram 
in Fig[9l In this case gauge invariance is not present (36| 
and reentrance do not occur, as expected. On the other 
side we still observe ferromagnetic and paramagnetic or- 
der, divided by the phase boundary (red thick line in Fig. 
[8]), where a multicritical point occur. As in the Q — 2,4 
case, any initial condition close to the phase boundary 
and above the multicritical point flows to an unstable 
finite temperature fix point, while any initial condition 
close to the boundary and below the multicritical point 
fiows to a strong coupling zero temperature distribution, 
even though gauge symmetry is not present. This implies 
that at small, finite values of the clock state parameter 
the strong violation of universality, discussed in the past 
for the (RBIM) \^ , also occurs in the clock state model, 
for values of the clock state parameter Q < 8. 

In the case of larger values of the clock state parame- 
ter, e.g. Q = 16, and still at 13 = 0, we do not observe 
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FIG. 9: The Phase Digram for the clock state parameter value 
Q = 8, without (DM) interactions and in the presence of bi- 
modal exchange interactions, with an antiferromagnetic bond 
concentration p. Gauge invariance is not present in this prob- 
lem, since (DM) interactions are set to zero. Reentrance is 
not expected in this case and the Nishimori symmetry line is 
not present. We observe however the presence of a multicrit- 
ical point. Any initial condition above the multicritical point 
and close to the phase boundary, flows to a finite tempera- 
ture unstable fix point, as seen in the inset above, while any 
initial condition below the multicritical point (indicated my 
the cross) and close to the boundary flows to a strong cou- 
pling, low temperature fix distribution, as in the second inset. 
The insets plot, as a function of the KG iteration number the 

Average Potential, defined as y^J^^ | ^(Q, qW- 



the recursion relations ^ to flow to a ferromagnetic fix 
point anymore, but rather quasi-fixed line behavior of the 
(BKT) type appears at finite values of the antiferromag- 
netic bond concentration, as seen in Fig. 1101 

We are considering in fact the discrete version of the 
XY spin glass problem. No gauge invariance nor reen- 
trance is expected in this case, and we observe effective 
algebraic order at low values of the bimodal exchange in- 
teraction, while a paramagnetic phase occurs at higher 
values of disorder and temperature. Clearly no spin-glass 
order is observed in two dimension everywhere in the 
phase diagram. 

Global phase diagrams for the XYSG and XYDM 
problems are being evaluated for values of the clock state 
parameter Q = 16, 32, 64 and will be presented elsewhere. 

We did not consider yet the final case where both 
randomness in the exchange and (DM) interactions are 
present, even though we are aware that this was the main 
question in order to discuss oscillator arrays, and leave 
this question for future work. 

We believe we have set up a consistent machinery to 
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FIG. 10: Number of iterations required to flow to the param- 
agnetic sink, as a function of temperature for the clock state 
parameter value Q = 16 for different values of the antiferro- 
magnetic bond concentration p = 0.0, p = 0.01, • ■ • , 0.07. The 
figure indicates that effective algebraic order is observed for 
finite, small values of the disorder strength in the exchange 
interaction. From the location of these lines it is possible 
to reconstruct the phase diagram of the XYSG problem at 
Q = 16. 

evaluate the phase diagram in the general case of large 
values of Q and disorder of both types and results in this 
direction wiU be able to answer on the issue of reentrance 
for continuous spin models, at least within the MK ap- 
proximation. 

When both disorder in the exchange and (DM) interac- 
tions is present, the renormalization group recursions in- 
volve a two dimensional probability distribution, so that 
the binning technique required is slightly more compli- 
cated than the one considered in the absence of (DM) 
interactions. A similar type of situation has been con- 
sidered for the Random Field Random Bond (RFRB) 
model, where three dimensional probability distributions 
were considered, and it is just a technical issue (but 
rather tedious) to extend the techniques developed in 
that case to the case of random exchange and random 
(DM) interactions, at finite, large values of the clock state 
parameter Q. 



VIII. CONCLUSIONS AND FUTURE 
PERSPECTIVES 

We conclude with a few considerations and predictions 
on what we expect to see in the general case of large 
values of the clock state parameter Q and disorder of the 
(DM) and exchange type, according to the preliminary 




Disorder Strength, A 



FIG. 11: Qualitative behaviour of the phase diagram in the 
case of the two dimensional ferromagnet with (DM) random 
interactions of width A. (A) Real space renormalization [2ll ]. 
suggesting reentrance that extends all the way down to van- 
ishing disorder strength. (B) Absence of reentrance and ver- 
tical boundary at finite values of the disorder strength as sug- 
gested, e.g. in reference [S^. (C) Reentrant behavior of the 
phase boundary and a (BKT) phase that extend at finite val- 
ues of the disorder strength for low temperature values. 



results we obtained so far. 

If both randomness in the exchange and (DM) interac- 
tion are considered in a way such that gauge invariance 
is recovered, we expect to observe a (BKT) phase and 
reentrant behavior of the phase boundary of the type 
discussed above. Note that the exchange randomness is 
usually considered as irrelevant, when (DM) random in- 
teractions are present, so similar conclusions should ap- 
ply to the (XYDM) problem. 

We expect the phase boundary to reenter below the 
Nishimori line, and finally to intercept the zero temper- 
ature axis with vertical slope, at finite, non — vanishing 
values of the disorder strength. This would imply that 
(AO) is stable against the intrinsic frequency distribution 
and it is reasonable to ask whether one can observe order 
of the (BKT) in arrays of two dimensional, non-identical 
oscillators, where technically one should consider tem- 
perature to be small. Most likely, the reason why pre- 
vious calculations within the (MK) approximation have 
not been able to show the reentrance predicted is simply 
because the recursion relations considered in [H, |3l| were 
not taking properly into account the (DM) interactions, 
as stressed above. This is likely the reason why the (MK) 
phase diagram discussed in [23] is not reentrant, as we 
also observe for the case of Q = 8 in Fig. pT|) where 
the (DM) interaction is explicitly set to zero and gauge 
invariance is not present. A careful analysis of the recur- 
sion relations ([6]) in the presence of (DM) interactions, 
and the corresponding phase diagram, both at small and 
large values of the clock state parameter Q is in progress 
and will be presented shortly. 
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The results we expect can be summarised as follows. 
The phase diagram is reentrant as predicted long ago, 
but reentrance (see curve B in Fig. [TT|) reduces and the 
phase boundary between the (BKT) phase and the para- 
magnetic phase intercept, with vertical slope, the zero 
temperature axis at finite values of the disorder strength. 
This intermediate scenario would be consistent both with 
the reentrant behavior predicted by classical renormaliza- 
tion group arguments [21[ and with the results indicating 
that the phase boundary intercepts the zero temperature 
axis with vertical slope [s^. 
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